########################################################################################################################################################
# September 24 2020
# Replication File
# Rebecca Cordell
# The Political Costs of Abusing Human Rights: International Cooperation in Extraordinary Rendition
# Journal of Conflict Resolution
# https://journals.sagepub.com/home/jcr  
########################################################################################################################################################

# Clear work environment
rm(list=ls())

# Required packages
#install.packages("arm")
#install.packages("stargazer")
#install.packages("lmtest")
#install.packages("pscl")
#install.packages("cem")
#install.packages("ggplot2")
#install.packages("reshape2")
#install.packages("ggpubr")
#install.packages("margins")
library(arm)
library(stargazer)
library(lmtest)
library(pscl)
library(cem)
library(ggplot2)
library(reshape2)
library(ggpubr)
library(margins)
options(scipen=999)

###############
# Read in data
###############

# Pre-matched data
rdi_data<-read.csv("Cordell_2020_PoliticalCosts.csv", header=TRUE)

# Post-matched data
matched_data<-rdi_data[rdi_data$matched==1,]

###############################################################
# Table 1: Probit Regression, Electoral Defeat, Matched Sample
###############################################################

# Model 1 Left Orientation Interaction
tab1_m1<- glm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= matched_data)

# Model 2 Interaction Components and Control Variables
tab1_m2<- glm(electoral.defeat ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Model 3 Full Model Left Orientation Interaction
tab1_m3<- glm(electoral.defeat ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(tab1_m1, tab1_m2, tab1_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(tab1_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(tab1_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(tab1_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(tab1_m1), digits=3) # McFadden
round(pR2(tab1_m2), digits=3) # McFadden
round(pR2(tab1_m3), digits=3) # McFadden

#####################################################################################
# Figure 1 and 2: Timeline of Revelation of Participation in RDI and Elections
#####################################################################################

# Structure data
df <- melt(matched_data, id.vars = "country", measure.vars = c("election.year", "revelation.year", "revelation.election.year"))
colnames(matched_data)
fig1_2<-matched_data[c(3,25)]
df<-merge(df, fig1_2, by=c("country"), all.x=TRUE)
df$variable <- factor(df$variable, levels=c("revelation.year", "election.year", "revelation.election.year"), labels=c("Revelation", "Election", "Revelation and Election"))
df$matched <- factor(df$matched)
df_1<-df[1:219,]
df_2<-df[220:435,]
fig1<-rdi_data[1:73,]
fig2<-rdi_data[74:145,]

# Figure 1 data
df_1 <- melt(fig1, id.vars = "country", measure.vars = c("election.year", "revelation.year", "revelation.election.year"))
fig1<-fig1[c(3,25)]
df_1<-merge(df_1, fig1, by=c("country"), all.x=TRUE)
df_1$variable <- factor(df_1$variable, levels=c("revelation.year", "election.year", "revelation.election.year"), labels=c("Revelation", "Election", "Revelation and Election"))
df_1$matched <- factor(df_1$matched)
df_1$value <- factor(df_1$value)
df_1$country <- factor(df_1$country,levels=rev(unique(df_1$country)))

# Figure 2 data
df_2 <- melt(fig2, id.vars = "country", measure.vars = c("election.year", "revelation.year", "revelation.election.year"))
fig2<-fig2[c(3,25)]
df_2<-merge(df_2, fig2, by=c("country"), all.x=TRUE)
df_2$variable <- factor(df_2$variable, levels=c("revelation.year", "election.year", "revelation.election.year"), labels=c("Revelation", "Election", "Revelation and Election"))
df_2$matched <- factor(df_2$matched)
df_2$value <- factor(df_2$value)
df_2$country <- factor(df_2$country,levels=rev(unique(df_2$country)))

# Figure 1: Timeline of Revelation of Participation in RDI and Elections (cont’d on next page)
fig1<-ggplot(data=subset(df_1, !is.na(value)), aes(y = country, x = value)) + labs(y="Country", x="Year")+
  geom_point(aes(group = variable, shape = variable, fill=matched)) +
  scale_shape_manual(values=c(24,21,22),labels = c("Revelation","Election", "Revelation and Election"))
fig1<-fig1+ labs(fill = "Matched Sample",shape = "Event")
fig1<-fig1+ guides(shape = guide_legend(reverse=TRUE))
fig1<-fig1+theme(axis.text.y=element_text(size=2))
fig1<-fig1+theme_bw() + theme(panel.border = element_blank(),
                              panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.text.y = element_text(size=6))
fig1<-fig1+scale_fill_manual(values=c("1"="black","0"="grey"),labels = c("Yes", "No"))
fig1<-fig1+guides(shape = guide_legend(order = 1), fill=guide_legend(order=2, override.aes=list(colour=c("1"="black","0"="grey"))))
fig1

# Figure 2: Timeline of Revelation of Participation in RDI and Elections (cont’d from previous page)
fig2<-ggplot(data=subset(df_2, !is.na(value)), aes(y = country, x = value)) + labs(y="Country", x="Year")+
  geom_point(aes(group = variable, shape = variable, fill=matched)) +
  scale_shape_manual(values=c(24,21,22),labels = c("Revelation","Election", "Revelation and Election"))
fig2<-fig2+ labs(fill = "Matched Sample",shape = "Event")
fig2<-fig2+ guides(shape = guide_legend(reverse=TRUE))
fig2<-fig2+theme(axis.text.y=element_text(size=2))
fig2<-fig2+theme_bw() + theme(panel.border = element_blank(),
                              panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.text.y = element_text(size=6))
fig2<-fig2+scale_fill_manual(values=c("1"="black","0"="grey"),labels = c("Yes", "No"))
fig2<-fig2+guides(shape = guide_legend(order = 1), fill=guide_legend(order=2, override.aes=list(colour=c("1"="black","0"="grey"))))
fig2

###################################################################################################
# Figure 3: Predicted Probabilities with 95% Confidence Intervals, Electoral Defeat, Matched Sample
###################################################################################################

# Model 3 Full Model Left Orientation Interaction
tab1_m3.pe <- summary(tab1_m3)$coefficients[,1]
tab1_m3.vcov <- vcov(tab1_m3)
Coefdraws.tab1_m3 <- MASS::mvrnorm(n = 10000, mu =  tab1_m3.pe, Sigma = tab1_m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r0l0 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r0l1 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r1l0 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r1l1 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mar = c(5,5,5,5))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c("Revelation & Left", "No Revelation & Left",
                                       "Revelation & Not Left", "No Revelation & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)

########################################################################################################################
# Figure 4: Bayesian Logistic Regression Coefficients, Electoral Defeat, Matched Sample, Democracies and Non-Democracies
########################################################################################################################

# Bayesian Logistic Regression, Electoral Defeat, Matched Sample, Democracies
democ<-matched_data[matched_data$democ==1,]

# Model 1 Left Orientation Interaction
taba5_m1<- bayesglm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= democ, prior.scale=2.5, prior.scale.for.intercept=10)

# Model 2 Interaction Components and Control Variables
taba5_m2<- bayesglm(electoral.defeat ~ revelation + left.orientation+  election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="logit"), data= democ, prior.scale=2.5, prior.scale.for.intercept=10)

# Model 3 Full Model Left Orientation Interaction
taba5_m3<- bayesglm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation  + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="logit"), prior.scale=2.5, prior.scale.for.intercept=10, data=democ)

# Bayesian Logistic Regression, Electoral Defeat, Matched Sample, Non-democracies
nondemoc<-matched_data[matched_data$nondemoc==1,]

# Model 1 Left Orientation Interaction
taba6_m1<- bayesglm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= nondemoc, prior.scale=2.5, prior.scale.for.intercept=10)

# Model 2 Interaction Components and Control Variables
taba6_m2<- bayesglm(electoral.defeat ~ revelation + left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="logit"), data= nondemoc, prior.scale=2.5, prior.scale.for.intercept=10)

# Model 3 Full Model Left Orientation Interaction
taba6_m3<- bayesglm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="logit"), prior.scale=2.5, prior.scale.for.intercept=10, data=nondemoc)

# Model 3 Full Model Left Orientation Interaction, Democracies
democ_slopes <- coef(summary(taba5_m3))[c("revelation", "left.orientation", "revelation:left.orientation"), 1]  #' democ_slopes
democ_ses <- coef(summary(taba5_m3))[c("revelation", "left.orientation", "revelation:left.orientation"), 2]  #' SEs
democ_coef<-cbind(democ_slopes,democ_ses)
democ_coef<-as.data.frame(democ_coef)
democ_ci_low<-democ_coef$democ_slopes - (qnorm(0.95) * democ_coef$democ_ses)
democ_ci_high<-democ_coef$democ_slopes + (qnorm(0.95) * democ_coef$democ_ses)
democ_ci<-cbind(democ_ci_low,democ_ci_high)
democ_coef<-cbind(democ_coef,democ_ci)
democ_coef$full.name<-c("Revelation", "Left Orientation", "Revelation*Left Orientation")
democ_coef$order<-c("c","b","a")
democ_coef$full.name <- factor(democ_coef$full.name, levels = democ_coef$full.name[order(democ_coef$order)])

# Democracies Figure
fig4_democ<-ggplot(data = democ_coef) +
  geom_point(aes(full.name, democ_slopes)) +
  geom_errorbar(aes(x = full.name, ymin = democ_ci_low, ymax = democ_ci_high), size=.3,width=.2) +
  xlab('Variable') +
  ylab('Estimate') +
  coord_flip() +
  geom_hline(yintercept = 0) +
  theme_classic(base_size = 10) +
  ggtitle("Democracies") +
  theme(plot.title=element_text(hjust=0.5), plot.margin = margin(10, 10, 10, 10)) +
  theme(axis.text.y = element_text(angle = 360)) +
  theme(axis.text.x = element_text(angle = 360)) 

# Model 3 Full Model Left Orientation Interaction, Non-democracies
nondemoc_slopes <- coef(summary(taba6_m3))[c("revelation", "left.orientation", "revelation:left.orientation"), 1]  #' nondemoc_slopes
nondemoc_ses <- coef(summary(taba6_m3))[c("revelation", "left.orientation", "revelation:left.orientation"), 2]  #' SEs
nondemoc_coef<-cbind(nondemoc_slopes,nondemoc_ses)
nondemoc_coef<-as.data.frame(nondemoc_coef)
nondemoc_ci_low<-nondemoc_coef$nondemoc_slopes - (qnorm(0.95) * nondemoc_coef$nondemoc_ses)
nondemoc_ci_high<-nondemoc_coef$nondemoc_slopes + (qnorm(0.95) * nondemoc_coef$nondemoc_ses)
nondemoc_ci<-cbind(nondemoc_ci_low,nondemoc_ci_high)
nondemoc_coef<-cbind(nondemoc_coef,nondemoc_ci)
nondemoc_coef$full.name<-c("Revelation", "Left Orientation", "Revelation*Left Orientation")
nondemoc_coef$order<-c("c","b","a")
nondemoc_coef$full.name <- factor(nondemoc_coef$full.name, levels = nondemoc_coef$full.name[order(nondemoc_coef$order)])

# Non-democracies Figure
fig4_nondemoc<-ggplot(data = nondemoc_coef) +
  geom_point(aes(full.name, nondemoc_slopes)) +
  geom_errorbar(aes(x = full.name, ymin = nondemoc_ci_low, ymax = nondemoc_ci_high), size=.3,width=.2) +
  xlab('Variable') +
  ylab('Estimate') +
  coord_flip() +
  geom_hline(yintercept = 0) +
  theme_classic(base_size = 10) +
  ggtitle("Non-democracies") +
  theme(plot.title=element_text(hjust=0.5), plot.margin = margin(10, 10, 10, 10)) +
  theme(axis.text.y = element_text(angle = 360)) +
  theme(axis.text.x = element_text(angle = 360)) 

# Create figure
ggarrange(fig4_democ,fig4_nondemoc, ncol=2, nrow=1)

#################################################################
# Results - In Text Quotations
#################################################################

# Specifically, the average marginal effect of Revelation on Electoral Defeat is 0.106 when Left Orientation is 0 (0.111 standard error) and 0.360 when it is 1 (0.099 standard error).
tab1_m3_me<-margins(tab1_m3, at = list(left.orientation = 0:1))
tab1_m3_me<-summary(tab1_m3_me)
tab1_m3_me<-tab1_m3_me[11:12,]

tab1_m3_me

# As expected, left of center governments that were caught participating in RDI were most likely to be removed from office, with a predicted probability of Electoral Defeat at 82%.
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))
ystar.r1l1 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)

round(yhat.sim.pe.r1l1, digits=2)*100

# In comparison, the predicted probability of Electoral Defeat for a left of center governments that was not caught participating in RDI is only 21%.
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))
ystar.r0l1 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)

round(yhat.sim.pe.r0l1, digits=2)*100

# This finding is striking as left of center governments were 61 percentage points more likely to lose at the election if they were caught cooperating in RDI.

round(yhat.sim.pe.r1l1-yhat.sim.pe.r0l1, digits=2)*100

# On the other hand, non-left of center governments were only 10% more likely to experience electoral defeat if they were caught participating in RDI than those that did not participate (with a predicted probability of 5% for the former and 14% for the latter).
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))
ystar.r1l0 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)

X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))
ystar.r0l0 <- Coefdraws.tab1_m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)

round(yhat.sim.pe.r1l0-yhat.sim.pe.r0l0, digits=2)*100
round(yhat.sim.pe.r1l0, digits=2)*100
round(yhat.sim.pe.r0l0, digits=2)*100

# The effect of Revelation*Left Orientation on Electoral Defeat is positive and significant at the 90% confidence level for the democracy subsample (with a p value of 0.519).
taba5_m3_coef<-summary(taba5_m3)
taba5_m3_coef<-taba5_m3_coef$coefficients
taba5_m3_coef<-taba5_m3_coef[4,]

taba5_m3_coef

# For the democracy sub-sample, the average marginal effect of Revelation on Electoral Defeat is 0.331 when Left Orientation is 0 (0.244 standard error) and 0.658 when it is 1 (0.154 standard error).
taba5_m3_me<-margins(taba5_m3, at = list(left.orientation = 0:1))
taba5_m3_me<-summary(taba5_m3_me)
taba5_m3_me<-taba5_m3_me[9:10,]

taba5_m3_me

# For the non-democracy sub-sample, the average marginal effect of Revelation on Electoral Defeat is -0.033 when Left Orientation is 0 (0.111 standard error) to 0.131 when it is 1 (0.179 standard error).
taba6_m3_me<-margins(taba6_m3, at = list(left.orientation = 0:1))
taba6_m3_me<-summary(taba6_m3_me)
taba6_m3_me<-taba6_m3_me[9:10,]

taba6_m3_me

# Specifically, the difference in predicted probabilities of Electoral Defeat for left of center governments that participated in RDI is 46% higher in democracies than in non-democracies (with a predicted probability of 89% for the former and 43% for the latter.
taba5_m3.pe <- summary(taba5_m3)$coefficients[,1]
taba5_m3.vcov <- vcov(taba5_m3)
Coefdraws.taba5_m3 <- MASS::mvrnorm(n = 10000, mu =  taba5_m3.pe, Sigma = taba5_m3.vcov)
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(democ$election.proximity),
                casualties.log= mean(democ$casualties.log),
                terrorism.log= mean(democ$terrorism.log),
                economic.crisis= mean(democ$economic.crisis),
                term.limit= mean(democ$term.limit))
ystar.r1l1 <- Coefdraws.taba5_m3 %*% t(as.matrix(X.r1l1))
democ.yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
democ.yhat.sim.pe.r1l1 <- apply(democ.yhat.sim.r1l1, 2, quantile, probs = 0.50)

taba6_m3.pe <- summary(taba6_m3)$coefficients[,1]
taba6_m3.vcov <- vcov(taba6_m3)
Coefdraws.taba6_m3 <- MASS::mvrnorm(n = 10000, mu =  taba6_m3.pe, Sigma = taba6_m3.vcov)
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(nondemoc$election.proximity),
                casualties.log= mean(nondemoc$casualties.log),
                terrorism.log= mean(nondemoc$terrorism.log),
                economic.crisis= mean(nondemoc$economic.crisis),
                term.limit= mean(nondemoc$term.limit))
ystar.r1l1 <- Coefdraws.taba6_m3 %*% t(as.matrix(X.r1l1))
nondemoc.yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
nondemoc.yhat.sim.pe.r1l1 <- apply(nondemoc.yhat.sim.r1l1, 2, quantile, probs = 0.50)

round(democ.yhat.sim.pe.r1l1-nondemoc.yhat.sim.pe.r1l1, digits=2)*100
round(democ.yhat.sim.pe.r1l1, digits=2)*100
round(nondemoc.yhat.sim.pe.r1l1, digits=2)*100

########################################################################################################################################################
# Appendices
########################################################################################################################################################

#########################################################################################
# Table A1: Countries that participated in post-9/11 RDI (Open Society Foundations, 2013)
#########################################################################################

# See Open Society Foundations (2013). Globalizing Torture: CIA Secret Detention and Extraordinary Rendition. Available at https://www.opensocietyfoundations.org/reports/globalizing-torture-cia-secret-detention-and-extraordinary-rendition

#######################################################################################################
# Table A2: Descriptive Statistics of Independent and Control Variables, after Coarsened Exact Matching
#######################################################################################################

# N
revelation_n<-length(matched_data$revelation)
left.orientation_n<-length(matched_data$left.orientation)
regime.type_n<-length(matched_data$regime.type)
election.proximity_n<-length(matched_data$election.proximity) 
casualties.log_n<-length(matched_data$casualties.log)  
terrorism.log_n<-length(matched_data$terrorism.log) 
economic.crisis_n<-length(matched_data$economic.crisis)  
term.limit_n<-length(matched_data$term.limit)

# Min, Mean, Max
revelation_min_mean_max<-summary(matched_data$revelation)
left.orientation_min_mean_max<-summary(matched_data$left.orientation)   
regime.type_min_mean_max<-summary(matched_data$regime.type)  
election.proximity_min_mean_max<-summary(matched_data$election.proximity) 
casualties.log_min_mean_max<-summary(matched_data$casualties.log)  
terrorism.log_min_mean_max<-summary(matched_data$terrorism.log) 
economic.crisis_min_mean_max<-summary(matched_data$economic.crisis)  
term.limit_min_mean_max<-summary(matched_data$term.limit)
revelation_mean<-revelation_min_mean_max[c(4)]
left.orientation_mean<-left.orientation_min_mean_max[c(4)] 
regime.type_mean<-regime.type_min_mean_max[c(4)] 
election.proximity_mean<-election.proximity_min_mean_max[c(4)]
casualties.log_mean<-casualties.log_min_mean_max[c(4)]
terrorism.log_mean<-terrorism.log_min_mean_max[c(4)]
economic.crisis_mean<-economic.crisis_min_mean_max[c(4)]
term.limit_mean<-term.limit_min_mean_max[c(4)]
revelation_min<-revelation_min_mean_max[c(1)]
left.orientation_min<-left.orientation_min_mean_max[c(1)] 
regime.type_min<-regime.type_min_mean_max[c(1)] 
election.proximity_min<-election.proximity_min_mean_max[c(1)]
casualties.log_min<-casualties.log_min_mean_max[c(1)]
terrorism.log_min<-terrorism.log_min_mean_max[c(1)]
economic.crisis_min<-economic.crisis_min_mean_max[c(1)]
term.limit_min<-term.limit_min_mean_max[c(1)]
revelation_max<-revelation_min_mean_max[c(6)]
left.orientation_max<-left.orientation_min_mean_max[c(6)] 
regime.type_max<-regime.type_min_mean_max[c(6)] 
election.proximity_max<-election.proximity_min_mean_max[c(6)]
casualties.log_max<-casualties.log_min_mean_max[c(6)]
terrorism.log_max<-terrorism.log_min_mean_max[c(6)]
economic.crisis_max<-economic.crisis_min_mean_max[c(6)]
term.limit_max<-term.limit_min_mean_max[c(6)]

# Standard Deviation
revelation_sd<-sd(matched_data$revelation)
left.orientation_sd<-sd(matched_data$left.orientation)   
regime.type_sd<-sd(matched_data$regime.type)  
election.proximity_sd<-sd(matched_data$election.proximity) 
casualties.log_sd<-sd(matched_data$casualties.log)  
terrorism.log_sd<-sd(matched_data$terrorism.log) 
economic.crisis_sd<-sd(matched_data$economic.crisis)  
term.limit_sd<-sd(matched_data$term.limit)

# Join Data
revelation_all<-cbind(revelation_n,revelation_mean,revelation_sd,revelation_min,revelation_max)
left.orientation_all<-cbind(left.orientation_n,left.orientation_mean,left.orientation_sd,left.orientation_min,left.orientation_max) 
regime.type_all<-cbind(regime.type_n,regime.type_mean,regime.type_sd,regime.type_min,regime.type_max) 
election.proximity_all<-cbind(election.proximity_n,election.proximity_mean,election.proximity_sd,election.proximity_min,election.proximity_max)
casualties.log_all<-cbind(casualties.log_n,casualties.log_mean,casualties.log_sd,casualties.log_min,casualties.log_max)
terrorism.log_all<-cbind(terrorism.log_n,terrorism.log_mean,terrorism.log_sd,terrorism.log_min,terrorism.log_max)
economic.crisis_all<-cbind(economic.crisis_n,economic.crisis_mean,economic.crisis_sd,economic.crisis_min,economic.crisis_max)
term.limit_all<-cbind(term.limit_n,term.limit_mean,term.limit_sd,term.limit_min,term.limit_max)

# Create Table
taba2<-round(rbind(revelation_all,left.orientation_all,regime.type_all,election.proximity_all,casualties.log_all,terrorism.log_all,economic.crisis_all,term.limit_all), digits=3)
taba2

########################################################################################
# Table A3: Matched Sample, Balanced on all Control Variables Included in the Full Model
########################################################################################

# Structure treatment and control group data
treatment<-matched_data[matched_data$revelation==1,]
control<-matched_data[matched_data$revelation==0,]
treat_name<-as.data.frame(treatment$country)
control_name<-control$country
treat_na<-as.data.frame(c(1:43))
colnames(treat_na)[colnames(treat_na)=="c(1:43)"] <- "country"
colnames(treat_name)[colnames(treat_name)=="treatment$country"] <- "country"
treat_na$country<-gsub("\\d","",treat_na$country)
treat_name<-rbind(treat_name,treat_na)

# Create Table
taba3<-cbind(treat_name,control_name)
names(taba3)[names(taba3) == "country"] <- "treatment.group"
names(taba3)[names(taba3) == "control_name"] <- "control.group"
taba3

####################################################################################################################
# Table A4: Percent Balance Improvement of Difference in Means for Covariates Across the Treatment and Control Group
####################################################################################################################

# Prematched data
prematch<-rdi_data[,c("revelation","left.orientation","regime.type","election.proximity","casualties.log","terrorism.log","economic.crisis","term.limit")]

# Imbalance statistics, prematched data
pre.imbalance <- imbalance(group=prematch$revelation,data=prematch,drop=c("revelation"))

# Postmatched data
postmatch<-matched_data[,c("revelation","left.orientation","regime.type","election.proximity","casualties.log","terrorism.log","economic.crisis","term.limit")]

# Imbalance statistics, postmatched data
post.imbalance <- imbalance(group=postmatch$revelation,data=postmatch,drop=c("revelation"))

# Create table
pre.diff.means<-abs(pre.imbalance$tab$statistic)
post.diff.means<-abs(post.imbalance$tab$statistic)
taba4<-cbind(pre.diff.means,post.diff.means)
taba4<-as.data.frame(taba4)
diff<-taba4$post.diff.means-taba4$pre.diff.means
taba4$percent.change<-abs(diff/taba4$pre.diff.means*100)
taba4<-round(taba4,digits=3)
taba4$change<-taba4$pre.diff.means-taba4$post.diff.means
taba4$change<-ifelse(taba4$change <0, "Increase", "Decrease")
taba4<-taba4[c(1,2,4,3)]
rownames(taba4)<-c("Left Orientation","Regime Type","Election Proximity","Casualties (log)","Terrorism (log)","Economic Crisis","Term Limit")
taba4

#######################################################################################
# Table A5: Bayesian Logistic Regression, Electoral Defeat, Matched Sample, Democracies
#######################################################################################

# See Line 224-246

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
taba5_m1_dta<-taba5_m1
taba5_m2_dta<-taba5_m2
taba5_m3_dta<-taba5_m3
class(taba5_m1_dta) <- c("glm","lm")
taba5_m1_dta$call[1] <- quote(lm())
class(taba5_m2_dta) <- c("glm","lm")
taba5_m2_dta$call[1] <- quote(lm())
class(taba5_m3_dta) <- c("glm","lm")
taba5_m3_dta$call[1] <- quote(lm())
stargazer(taba5_m1_dta, taba5_m2_dta, taba5_m3_dta, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba5_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba5_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba5_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba5_m1), digits=3) # McFadden
round(pR2(taba5_m2), digits=3) # McFadden
round(pR2(taba5_m3), digits=3) # McFadden

###########################################################################################
# Table A6: Bayesian Logistic Regression, Electoral Defeat, Matched Sample, Non-democracies
###########################################################################################

# See Line 215-237

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
taba6_m1_dta<-taba6_m1
taba6_m2_dta<-taba6_m2
taba6_m3_dta<-taba6_m3
class(taba6_m1_dta) <- c("glm","lm")
taba6_m1_dta$call[1] <- quote(lm())
class(taba6_m2_dta) <- c("glm","lm")
taba6_m2_dta$call[1] <- quote(lm())
class(taba6_m3_dta) <- c("glm","lm")
taba6_m3_dta$call[1] <- quote(lm())
stargazer(taba6_m1_dta, taba6_m2_dta, taba6_m3_dta, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba6_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba6_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba6_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba6_m1), digits=3) # McFadden
round(pR2(taba6_m2), digits=3) # McFadden
round(pR2(taba6_m3), digits=3) # McFadden

##############################################################################
# Table A7: Probit Regression, Source of Leader Support Change, Matched Sample
##############################################################################

# Model 1 Left Orientation Interaction
taba7_m1<- glm(solschangedum_any_rev ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= matched_data)

# Model 2 Interaction Components and Control Variables
taba7_m2<- glm(solschangedum_any_rev ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Model 3 Full Model Left Orientation Interaction
taba7_m3<- glm(solschangedum_any_rev ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba7_m1, taba7_m2, taba7_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba7_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba7_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba7_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba7_m1), digits=3) # McFadden
round(pR2(taba7_m2), digits=3) # McFadden
round(pR2(taba7_m3), digits=3) # McFadden

###################################################################################################################################
# Table A8 and A9: Probit Regression, Electoral Defeat, Excluding Countries that Cooperated and were Left of Center, Matched Sample
###################################################################################################################################

# Prepare Data (excluding each country)
alg<-matched_data[matched_data$country!="Algeria" ,] 
can<-matched_data[matched_data$country!="Canada",]      
cz<-matched_data[matched_data$country!="Czech Republic" ,] 
ger<-matched_data[matched_data$country!="Germany" ,]
mac<-matched_data[matched_data$country!="Macedonia"   ,]    
pol<-matched_data[matched_data$country!="Poland"   ,]   
por<-matched_data[matched_data$country!="Portugal" ,]             
swe<-matched_data[matched_data$country!="Sweden"   ,]               
uzb<-matched_data[matched_data$country!="Uzbekistan" ,]
sou<-matched_data[matched_data$country!="South Africa"   ,]               
lib<-matched_data[matched_data$country!="Libya" ,]

# Table A8: Probit Regression, Electoral Defeat, Excluding Countries that Cooperated and were Left of Center, Matched Sample (cont’d on next page)

# Model 3 Full Model Left Orientation Interaction
taba8.alg.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= alg)
taba8.can.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= can)
taba8.cz.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= cz)
taba8.ger.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= ger)
taba8.lib.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= lib)
taba8.mac.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= mac)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba8.alg.m3, taba8.can.m3, taba8.cz.m3, taba8.ger.m3, taba8.lib.m3, taba8.mac.m3, type = "text", order=c(1,2,9,3,4,5,6,7,8))

# Table LR chi2 and Prob>chi2
round(lrtest(taba8.alg.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba8.can.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba8.cz.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba8.ger.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba8.lib.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba8.mac.m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba8.alg.m3), digits=3) # McFadden
round(pR2(taba8.can.m3), digits=3) # McFadden
round(pR2(taba8.cz.m3), digits=3) # McFadden
round(pR2(taba8.ger.m3), digits=3) # McFadden
round(pR2(taba8.lib.m3), digits=3) # McFadden
round(pR2(taba8.mac.m3), digits=3) # McFadden

# Table A9: Probit Regression, Electoral Defeat, Excluding Countries that Cooperated and were Left of Center, Matched Sample (cont’d from previous page)
taba9.pol.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= pol)
taba9.por.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= por)
taba9.sou.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= sou)
taba9.swe.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= swe)
taba9.uzb.m3<- glm(electoral.defeat ~ revelation + left.orientation +  revelation*left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= uzb)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
taba9<-stargazer(taba9.pol.m3, taba9.por.m3, taba9.sou.m3, taba9.swe.m3, taba9.uzb.m3, type = "text", order=c(1,2,9,3,4,5,6,7,8))

# Table LR chi2 and Prob>chi2
round(lrtest(taba9.pol.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba9.por.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba9.sou.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba9.swe.m3), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba9.uzb.m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba9.pol.m3), digits=3) # McFadden
round(pR2(taba9.por.m3), digits=3) # McFadden
round(pR2(taba9.sou.m3), digits=3) # McFadden
round(pR2(taba9.swe.m3), digits=3) # McFadden
round(pR2(taba9.uzb.m3), digits=3) # McFadden

#############################################################################################
# Table A10: Probit Regression, Electoral Defeat, Right of Center Governments, Matched Sample
#############################################################################################

# Model 1 Left Orientation Interaction
taba10_m1<- glm(electoral.defeat ~ revelation + right.orient + revelation*right.orient , family=binomial(link="probit"), data= matched_data)

# Model 2 Interaction Components and Control Variables
taba10_m2<- glm(electoral.defeat ~ revelation + right.orient + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Model 3 Full Model Left Orientation Interaction
taba10_m3<- glm(electoral.defeat ~ revelation + right.orient +  regime.type + revelation*right.orient + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba10_m1, taba10_m2, taba10_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba10_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba10_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba10_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba10_m1), digits=3) # McFadden
round(pR2(taba10_m2), digits=3) # McFadden
round(pR2(taba10_m3), digits=3) # McFadden

#################################################################################################################################################################
# Table A11: Probit Regression, Electoral Defeat, Excluding Countries with an Islamic Majority Population that Cooperated and were Left of Center, Matched Sample
#################################################################################################################################################################

# Exclude US, Algeria and Libya
no_us<-matched_data[matched_data$country!="United States of America",]
islam<-no_us[!(no_us$country == "Algeria" | no_us$country == "Libya"),]

# Model 1 Left Orientation Interaction
taba11_m1<- glm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), islam)

# Model 2 Interaction Components and Control Variables
taba11_m2<- glm(electoral.defeat ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), islam)

# Model 3 Full Model Left Orientation Interaction
taba11_m3<- glm(electoral.defeat ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), islam)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba11_m1, taba11_m2, taba11_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba11_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba11_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba11_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba11_m1), digits=3) # McFadden
round(pR2(taba11_m2), digits=3) # McFadden
round(pR2(taba11_m3), digits=3) # McFadden

##############################################################################################################################################################################
# Table A12: Probit Regression, Electoral Defeat, Excluding Countries with above Average Number of U.S. Military Bases that Cooperated and were Left of Center, Matched Sample 
##############################################################################################################################################################################

# Calculate mean value for number of U.S. bases
mean(no_us$us.bases)
bases<-no_us[!(no_us$revelation == 1 & no_us$left.orientation == 1 & no_us$us.bases>=195.6894),]

# Model 1 Left Orientation Interaction
taba12_m1<- glm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= bases)

# Model 2 Interaction Components and Control Variables
taba12_m2<- glm(electoral.defeat ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= bases)

# Model 3 Full Model Left Orientation Interaction
taba12_m3<- glm(electoral.defeat ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= bases)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba12_m1, taba12_m2, taba12_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba12_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba12_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba12_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba12_m1), digits=3) # McFadden
round(pR2(taba12_m2), digits=3) # McFadden
round(pR2(taba12_m3), digits=3) # McFadden

###############################################################################################################################################################
# Table A13: Probit Regression, Electoral Defeat, Excluding Countries with above Average US Trade%GDPPC that Cooperated and were Left of Center, Matched Sample
###############################################################################################################################################################

# Calculate mean value for number of U.S. Trade%GDPPC
mean(no_us$us.trade.gdp, na.rm=TRUE)
trade<-no_us[!(no_us$revelation == 1 & no_us$left.orientation == 1 & no_us$us.trade.gdp>=1540584),]

# Model 1 Left Orientation Interaction
taba13_m1<- glm(electoral.defeat ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= trade)

# Model 2 Interaction Components and Control Variables
taba13_m2<- glm(electoral.defeat ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= trade)

# Model 3 Full Model Left Orientation Interaction
taba13_m3<- glm(electoral.defeat ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= trade)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba13_m1, taba13_m2, taba13_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba13_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba13_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba13_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba13_m1), digits=3) # McFadden
round(pR2(taba13_m2), digits=3) # McFadden
round(pR2(taba13_m3), digits=3) # McFadden

#####################################################################################
# Table A14: Probit Regression, Electoral Defeat, Election Circa 1996, Matched Sample
#####################################################################################

# Model 1 Left Orientation Interaction
taba14_m1<- glm(electoral.defeat.1996 ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= matched_data)

# Model 2 Interaction Components and Control Variables
taba14_m2<- glm(electoral.defeat.1996 ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Model 3 Full Model Left Orientation Interaction
taba14_m3<- glm(electoral.defeat.1996 ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= matched_data)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba14_m1, taba14_m2, taba14_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba14_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba14_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba14_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba14_m1), digits=3) # McFadden
round(pR2(taba14_m2), digits=3) # McFadden
round(pR2(taba14_m3), digits=3) # McFadden

################################################################################################
# Table A15: Probit Regression, Electoral Defeat, Election Before the Revelation, Matched Sample
################################################################################################

prior<-matched_data[!is.na(matched_data$electoral.defeat.prior),]

# Model 1 Left Orientation Interaction
taba15_m1<- glm(electoral.defeat.prior ~ revelation + left.orientation + revelation*left.orientation , family=binomial(link="probit"), data= prior)

# Model 2 Interaction Components and Control Variables
taba15_m2<- glm(electoral.defeat.prior ~ revelation + left.orientation + regime.type + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= prior)

# Model 3 Full Model Left Orientation Interaction
taba15_m3<- glm(electoral.defeat.prior ~ revelation + left.orientation +  regime.type + revelation*left.orientation + election.proximity + casualties.log + terrorism.log + economic.crisis + term.limit, family=binomial(link="probit"), data= prior)

# Create Table Coefficient Estimates, Standard Errors, Observations, Log Likelihood, Akaike Inf. Crit.
stargazer(taba15_m1, taba15_m2, taba15_m3, type = "text")

# Table LR chi2 and Prob>chi2
round(lrtest(taba15_m1), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba15_m2), digits=3) # Chisq and Pr(>Chisq)
round(lrtest(taba15_m3), digits=3) # Chisq and Pr(>Chisq)

# Table Pseudo R2
round(pR2(taba15_m1), digits=3) # McFadden
round(pR2(taba15_m2), digits=3) # McFadden
round(pR2(taba15_m3), digits=3) # McFadden

#############################################################################################################
# Figure A1: Predicted Probabilities with 95% Confidence Intervals, Source of Leader Support Change
#############################################################################################################

# Model 3 Full Model Left Orientation Interaction
taba7_m3.pe <- summary(taba7_m3)$coefficients[,1]
taba7_m3.vcov <- vcov(taba7_m3)
Coefdraws.taba7_m3 <- MASS::mvrnorm(n = 10000, mu =  taba7_m3.pe, Sigma = taba7_m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r0l0 <- Coefdraws.taba7_m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r0l1 <- Coefdraws.taba7_m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r1l0 <- Coefdraws.taba7_m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(matched_data$election.proximity),
                regime.type= mean(matched_data$regime.type),
                casualties.log= mean(matched_data$casualties.log),
                terrorism.log= mean(matched_data$terrorism.log),
                economic.crisis= mean(matched_data$economic.crisis),
                term.limit= mean(matched_data$term.limit))

ystar.r1l1 <- Coefdraws.taba7_m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mar = c(5,5,5,5))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c("Revelation & Left", "No Revelation & Left",
                                       "Revelation & Not Left", "No Revelation & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)

#################################################################################################################################################################################################
# Figure A2: Predicted Probabilities with 95% Confidence Intervals, Electoral Defeat, Excluding Countries that Cooperated and were Left of Center, Matched Sample (cont’d on next page)
#################################################################################################################################################################################################

# Excluding Algeria 

# Model 3 Full Model Left Orientation Interaction
taba8.alg.m3.pe <- summary(taba8.alg.m3)$coefficients[,1]
taba8.alg.m3.vcov <- vcov(taba8.alg.m3)
Coefdraws.taba8.alg.m3 <- MASS::mvrnorm(n = 10000, mu =  taba8.alg.m3.pe, Sigma = taba8.alg.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r0l0 <- Coefdraws.taba8.alg.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r0l1 <- Coefdraws.taba8.alg.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r1l0 <- Coefdraws.taba8.alg.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r1l1 <- Coefdraws.taba8.alg.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mfrow=c(3,2),mgp=c(5,2,0),mar=c(4.1, 4.1, 2.1, 2.1))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Algeria")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
alg.plot <- recordPlot()

# Excluding Canada

# Model 3 Full Model Left Orientation Interaction
taba8.can.m3.pe <- summary(taba8.can.m3)$coefficients[,1]
taba8.can.m3.vcov <- vcov(taba8.can.m3)
Coefdraws.taba8.can.m3 <- MASS::mvrnorm(n = 10000, mu =  taba8.can.m3.pe, Sigma = taba8.can.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r0l0 <- Coefdraws.taba8.can.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r0l1 <- Coefdraws.taba8.can.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r1l0 <- Coefdraws.taba8.can.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r1l1 <- Coefdraws.taba8.can.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Canada")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)

can.plot <- recordPlot()

# Excluding Czech Republic

# Model 3 Full Model Left Orientation Interaction
taba8.cz.m3.pe <- summary(taba8.cz.m3)$coefficients[,1]
taba8.cz.m3.vcov <- vcov(taba8.cz.m3)
Coefdraws.taba8.cz.m3 <- MASS::mvrnorm(n = 10000, mu =  taba8.cz.m3.pe, Sigma = taba8.cz.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r0l0 <- Coefdraws.taba8.cz.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r0l1 <- Coefdraws.taba8.cz.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r1l0 <- Coefdraws.taba8.cz.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r1l1 <- Coefdraws.taba8.cz.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Czech Republic")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
cz.plot <- recordPlot()

# Excluding Germany 

# Model 3 Full Model Left Orientation Interaction
taba8.ger.m3.pe <- summary(taba8.ger.m3)$coefficients[,1]
taba8.ger.m3.vcov <- vcov(taba8.ger.m3)
Coefdraws.taba8.ger.m3 <- MASS::mvrnorm(n = 10000, mu =  taba8.ger.m3.pe, Sigma = taba8.ger.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r0l0 <- Coefdraws.taba8.ger.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r0l1 <- Coefdraws.taba8.ger.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r1l0 <- Coefdraws.taba8.ger.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r1l1 <- Coefdraws.taba8.ger.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "", main="Excluding Germany")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
ger.plot <- recordPlot()

# Excluding Libya 

# Model 3 Full Model Left Orientation Interaction
taba8.lib.m3.pe <- summary(taba8.lib.m3)$coefficients[,1]
taba8.lib.m3.vcov <- vcov(taba8.lib.m3)
Coefdraws.taba8.lib.m3 <- MASS::mvrnorm(n = 10000, mu =  taba8.lib.m3.pe, Sigma = taba8.lib.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r0l0 <- Coefdraws.taba8.lib.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r0l1 <- Coefdraws.taba8.lib.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r1l0 <- Coefdraws.taba8.lib.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r1l1 <- Coefdraws.taba8.lib.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Libya")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
lib.plot <- recordPlot()

# Excluding Macedonia

# Model 3 Full Model Left Orientation Interaction
taba8.mac.m3.pe <- summary(taba8.mac.m3)$coefficients[,1]
taba8.mac.m3.vcov <- vcov(taba8.mac.m3)
Coefdraws.taba8.mac.m3 <- MASS::mvrnorm(n = 10000, mu =  taba8.mac.m3.pe, Sigma = taba8.mac.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(mac$election.proximity),
                regime.type= mean(mac$regime.type),
                casualties.log= mean(mac$casualties.log),
                terrorism.log= mean(mac$terrorism.log),
                economic.crisis= mean(mac$economic.crisis), term.limit= mean(mac$term.limit))

ystar.r0l0 <- Coefdraws.taba8.mac.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(mac$election.proximity),
                regime.type= mean(mac$regime.type),
                casualties.log= mean(mac$casualties.log),
                terrorism.log= mean(mac$terrorism.log),
                economic.crisis= mean(mac$economic.crisis), term.limit= mean(mac$term.limit))

ystar.r0l1 <- Coefdraws.taba8.mac.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(mac$election.proximity),
                regime.type= mean(mac$regime.type),
                casualties.log= mean(mac$casualties.log),
                terrorism.log= mean(mac$terrorism.log),
                economic.crisis= mean(mac$economic.crisis), term.limit= mean(mac$term.limit))

ystar.r1l0 <- Coefdraws.taba8.mac.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(mac$election.proximity),
                regime.type= mean(mac$regime.type),
                casualties.log= mean(mac$casualties.log),
                terrorism.log= mean(mac$terrorism.log),
                economic.crisis= mean(mac$economic.crisis), term.limit= mean(mac$term.limit))

ystar.r1l1 <- Coefdraws.taba8.mac.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Macedonia")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
mac.plot <- recordPlot()

#######################################################################################################################################################################################################
# Figure A3: Predicted Probabilities with 95% Confidence Intervals, Electoral Defeat, Excluding Countries that Cooperated and were Left of Center, Matched Sample (cont’d from previous page)
#######################################################################################################################################################################################################

# Excluding Poland

# Model 3 Full Model Left Orientation Interaction
taba9.pol.m3.pe <- summary(taba9.pol.m3)$coefficients[,1]
taba9.pol.m3.vcov <- vcov(taba9.pol.m3)
Coefdraws.taba9.pol.m3 <- MASS::mvrnorm(n = 10000, mu =  taba9.pol.m3.pe, Sigma = taba9.pol.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r0l0 <- Coefdraws.taba9.pol.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r0l1 <- Coefdraws.taba9.pol.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r1l0 <- Coefdraws.taba9.pol.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(alg$election.proximity),
                regime.type= mean(alg$regime.type),
                casualties.log= mean(alg$casualties.log),
                terrorism.log= mean(alg$terrorism.log),
                economic.crisis= mean(alg$economic.crisis), term.limit= mean(alg$term.limit))

ystar.r1l1 <- Coefdraws.taba9.pol.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mfrow=c(3,2),mgp=c(5,2,0),mar=c(4.1, 4.1, 2.1, 2.1))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Poland")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
pol.plot <- recordPlot()

# Excluding Portugal

# Model 3 Full Model Left Orientation Interaction
taba9.por.m3.pe <- summary(taba9.por.m3)$coefficients[,1]
taba9.por.m3.vcov <- vcov(taba9.por.m3)
Coefdraws.taba9.por.m3 <- MASS::mvrnorm(n = 10000, mu =  taba9.por.m3.pe, Sigma = taba9.por.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r0l0 <- Coefdraws.taba9.por.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r0l1 <- Coefdraws.taba9.por.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r1l0 <- Coefdraws.taba9.por.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(can$election.proximity),
                regime.type= mean(can$regime.type),
                casualties.log= mean(can$casualties.log),
                terrorism.log= mean(can$terrorism.log),
                economic.crisis= mean(can$economic.crisis), term.limit= mean(can$term.limit))

ystar.r1l1 <- Coefdraws.taba9.por.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Portugal")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)

por.plot <- recordPlot()

# Excluding South Africa

# Model 3 Full Model Left Orientation Interaction
taba9.sou.m3.pe <- summary(taba9.sou.m3)$coefficients[,1]
taba9.sou.m3.vcov <- vcov(taba9.sou.m3)
Coefdraws.taba9.sou.m3 <- MASS::mvrnorm(n = 10000, mu =  taba9.sou.m3.pe, Sigma = taba9.sou.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r0l0 <- Coefdraws.taba9.sou.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r0l1 <- Coefdraws.taba9.sou.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r1l0 <- Coefdraws.taba9.sou.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(cz$election.proximity),
                regime.type= mean(cz$regime.type),
                casualties.log= mean(cz$casualties.log),
                terrorism.log= mean(cz$terrorism.log),
                economic.crisis= mean(cz$economic.crisis), term.limit= mean(cz$term.limit))

ystar.r1l1 <- Coefdraws.taba9.sou.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding South Africa")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
sou.plot <- recordPlot()

# Excluding Sweden 

# Model 3 Full Model Left Orientation Interaction
taba9.swe.m3.pe <- summary(taba9.swe.m3)$coefficients[,1]
taba9.swe.m3.vcov <- vcov(taba9.swe.m3)
Coefdraws.taba9.swe.m3 <- MASS::mvrnorm(n = 10000, mu =  taba9.swe.m3.pe, Sigma = taba9.swe.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r0l0 <- Coefdraws.taba9.swe.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r0l1 <- Coefdraws.taba9.swe.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r1l0 <- Coefdraws.taba9.swe.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(ger$election.proximity),
                regime.type= mean(ger$regime.type),
                casualties.log= mean(ger$casualties.log),
                terrorism.log= mean(ger$terrorism.log),
                economic.crisis= mean(ger$economic.crisis), term.limit= mean(ger$term.limit))

ystar.r1l1 <- Coefdraws.taba9.swe.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "", main="Excluding Sweden")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
swe.plot <- recordPlot()

# Excluding Uzbekistan 

# Model 3 Full Model Left Orientation Interaction
taba9.uzb.m3.pe <- summary(taba9.uzb.m3)$coefficients[,1]
taba9.uzb.m3.vcov <- vcov(taba9.uzb.m3)
Coefdraws.taba9.uzb.m3 <- MASS::mvrnorm(n = 10000, mu =  taba9.uzb.m3.pe, Sigma = taba9.uzb.m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r0l0 <- Coefdraws.taba9.uzb.m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r0l1 <- Coefdraws.taba9.uzb.m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r1l0 <- Coefdraws.taba9.uzb.m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(lib$election.proximity),
                regime.type= mean(lib$regime.type),
                casualties.log= mean(lib$casualties.log),
                terrorism.log= mean(lib$terrorism.log),
                economic.crisis= mean(lib$economic.crisis), term.limit= mean(lib$term.limit))

ystar.r1l1 <- Coefdraws.taba9.uzb.m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "",main="Excluding Uzbekistan")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c(" Revelation \n  & Left", " No Revelation \n  & Left",
                                       " Revelation \n  & Not Left", " No Revelation \n  & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)
uzb.plot <- recordPlot()

###############################################################################################################################################################################################
# Figure A4: Predicted Probabilities with 95% Confidence Intervals, Electoral Defeat, Excluding Countries with an Islamic Majority Population that Cooperated and were Left of Center 
###############################################################################################################################################################################################

# Model 3 Full Model Left Orientation Interaction
taba11_m3.pe <- summary(taba11_m3)$coefficients[,1]
taba11_m3.vcov <- vcov(taba11_m3)
Coefdraws.taba11_m3 <- MASS::mvrnorm(n = 10000, mu =  taba11_m3.pe, Sigma = taba11_m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(islam$election.proximity),
                regime.type= mean(islam$regime.type),
                casualties.log= mean(islam$casualties.log),
                terrorism.log= mean(islam$terrorism.log),
                economic.crisis= mean(islam$economic.crisis),
                term.limit= mean(islam$term.limit))

ystar.r0l0 <- Coefdraws.taba11_m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(islam$election.proximity),
                regime.type= mean(islam$regime.type),
                casualties.log= mean(islam$casualties.log),
                terrorism.log= mean(islam$terrorism.log),
                economic.crisis= mean(islam$economic.crisis),
                term.limit= mean(islam$term.limit))

ystar.r0l1 <- Coefdraws.taba11_m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(islam$election.proximity),
                regime.type= mean(islam$regime.type),
                casualties.log= mean(islam$casualties.log),
                terrorism.log= mean(islam$terrorism.log),
                economic.crisis= mean(islam$economic.crisis),
                term.limit= mean(islam$term.limit))

ystar.r1l0 <- Coefdraws.taba11_m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(islam$election.proximity),
                regime.type= mean(islam$regime.type),
                casualties.log= mean(islam$casualties.log),
                terrorism.log= mean(islam$terrorism.log),
                economic.crisis= mean(islam$economic.crisis),
                term.limit= mean(islam$term.limit))

ystar.r1l1 <- Coefdraws.taba11_m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mfrow=c(1,1), mar = c(5,5,5,5))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c("Revelation & Left", "No Revelation & Left",
                                       "Revelation & Not Left", "No Revelation & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)

############################################################################################################################################################################################################
# Figure A5: Predicted Probabilities with 95% Confidence Intervals, Electoral Defeat, Excluding Countries with above Average Number of U.S. Military Bases that Cooperated and were Left of Center 
############################################################################################################################################################################################################

# Model 3 Full Model Left Orientation Interaction
taba12_m3.pe <- summary(taba12_m3)$coefficients[,1]
taba12_m3.vcov <- vcov(taba12_m3)
Coefdraws.taba12_m3 <- MASS::mvrnorm(n = 10000, mu =  taba12_m3.pe, Sigma = taba12_m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(bases$election.proximity),
                regime.type= mean(bases$regime.type),
                casualties.log= mean(bases$casualties.log),
                terrorism.log= mean(bases$terrorism.log),
                economic.crisis= mean(bases$economic.crisis),
                term.limit= mean(bases$term.limit))

ystar.r0l0 <- Coefdraws.taba12_m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(bases$election.proximity),
                regime.type= mean(bases$regime.type),
                casualties.log= mean(bases$casualties.log),
                terrorism.log= mean(bases$terrorism.log),
                economic.crisis= mean(bases$economic.crisis),
                term.limit= mean(bases$term.limit))

ystar.r0l1 <- Coefdraws.taba12_m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(bases$election.proximity),
                regime.type= mean(bases$regime.type),
                casualties.log= mean(bases$casualties.log),
                terrorism.log= mean(bases$terrorism.log),
                economic.crisis= mean(bases$economic.crisis),
                term.limit= mean(bases$term.limit))

ystar.r1l0 <- Coefdraws.taba12_m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(bases$election.proximity),
                regime.type= mean(bases$regime.type),
                casualties.log= mean(bases$casualties.log),
                terrorism.log= mean(bases$terrorism.log),
                economic.crisis= mean(bases$economic.crisis),
                term.limit= mean(bases$term.limit))

ystar.r1l1 <- Coefdraws.taba12_m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mar = c(5,5,5,5))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c("Revelation & Left", "No Revelation & Left",
                                       "Revelation & Not Left", "No Revelation & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)

#############################################################################################################################################################################################
# Figure A6: Predicted Probabilities with 95% Confidence Intervals, Electoral Defeat, Excluding Countries with above Average US Trade%GDPPC that Cooperated and were Left of Center
#############################################################################################################################################################################################

# Model 3 Full Model Left Orientation Interaction
taba13_m3.pe <- summary(taba13_m3)$coefficients[,1]
taba13_m3.vcov <- vcov(taba13_m3)
Coefdraws.taba13_m3 <- MASS::mvrnorm(n = 10000, mu =  taba13_m3.pe, Sigma = taba13_m3.vcov)

# No Revelation & Not Left
X.r0l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(trade$election.proximity),
                regime.type= mean(trade$regime.type),
                casualties.log= mean(trade$casualties.log),
                terrorism.log= mean(trade$terrorism.log),
                economic.crisis= mean(trade$economic.crisis),
                term.limit= mean(trade$term.limit))

ystar.r0l0 <- Coefdraws.taba13_m3 %*% t(as.matrix(X.r0l0))
yhat.sim.r0l0 <- exp(ystar.r0l0) / (exp(ystar.r0l0) + 1)
yhat.sim.pe.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l0 <- apply(yhat.sim.r0l0, 2, quantile, probs = 0.975)

# No Revelation & Left
X.r0l1 <- cbind(intercept = 1, new_interactr0l1 = 1, new_interactr1l0 = 0, new_interactr1l1 = 0, 
                election.proximity= mean(trade$election.proximity),
                regime.type= mean(trade$regime.type),
                casualties.log= mean(trade$casualties.log),
                terrorism.log= mean(trade$terrorism.log),
                economic.crisis= mean(trade$economic.crisis),
                term.limit= mean(trade$term.limit))

ystar.r0l1 <- Coefdraws.taba13_m3 %*% t(as.matrix(X.r0l1))
yhat.sim.r0l1 <- exp(ystar.r0l1) / (exp(ystar.r0l1) + 1)
yhat.sim.pe.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r0l1 <- apply(yhat.sim.r0l1, 2, quantile, probs = 0.975)

# Revelation & Not Left
X.r1l0 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 1, new_interactr1l1 = 0, 
                election.proximity= mean(trade$election.proximity),
                regime.type= mean(trade$regime.type),
                casualties.log= mean(trade$casualties.log),
                terrorism.log= mean(trade$terrorism.log),
                economic.crisis= mean(trade$economic.crisis),
                term.limit= mean(trade$term.limit))

ystar.r1l0 <- Coefdraws.taba13_m3 %*% t(as.matrix(X.r1l0))
yhat.sim.r1l0 <- exp(ystar.r1l0) / (exp(ystar.r1l0) + 1)
yhat.sim.pe.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l0 <- apply(yhat.sim.r1l0, 2, quantile, probs = 0.975)

# Revelation & Left
X.r1l1 <- cbind(intercept = 1, new_interactr0l1 = 0, new_interactr1l0 = 0, new_interactr1l1 = 1, 
                election.proximity= mean(trade$election.proximity),
                regime.type= mean(trade$regime.type),
                casualties.log= mean(trade$casualties.log),
                terrorism.log= mean(trade$terrorism.log),
                economic.crisis= mean(trade$economic.crisis),
                term.limit= mean(trade$term.limit))

ystar.r1l1 <- Coefdraws.taba13_m3 %*% t(as.matrix(X.r1l1))
yhat.sim.r1l1 <- exp(ystar.r1l1) / (exp(ystar.r1l1) + 1)
yhat.sim.pe.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.50)
yhat.sim.lb.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.025)
yhat.sim.ub.r1l1 <- apply(yhat.sim.r1l1, 2, quantile, probs = 0.975)

# Create Figure
par(mar = c(5,5,5,5))
x.ax <- c(1, 2, 3, 4)
y.ax <- c(yhat.sim.pe.r1l1, yhat.sim.pe.r0l1,yhat.sim.pe.r1l0, yhat.sim.pe.r0l0)
plot(x.ax, y.ax, ylim = c(0,1), pch = 16, axes = F, ylab = "Predicted Probabilities", xlab = "")
axis(2, at = c(0, 0.25, 0.50, 0.75 , 1))
axis(1, at = c(1, 2, 3, 4), labels = c("Revelation & Left", "No Revelation & Left",
                                       "Revelation & Not Left", "No Revelation & Not Left"))
segments(x0 = 1, x1 = 1, y0 = yhat.sim.lb.r1l1, y1 = yhat.sim.ub.r1l1)
text(x = 1 , y = yhat.sim.lb.r1l1, labels = "-", cex = 1.5)
text(x = 1 , y = yhat.sim.ub.r1l1, labels = "-", cex = 1.5)
segments(x0 = 2, x1 = 2, y0 = yhat.sim.lb.r0l1, y1 = yhat.sim.ub.r0l1 )
text(x = 2 , y = yhat.sim.lb.r0l1, labels = "-", cex = 1.5)
text(x = 2 , y = yhat.sim.ub.r0l1, labels = "-", cex = 1.5)
segments(x0 = 3, x1 = 3, y0 = yhat.sim.lb.r1l0, y1 = yhat.sim.ub.r1l0 )
text(x = 3 , y = yhat.sim.lb.r1l0, labels = "-", cex = 1.5)
text(x = 3 , y = yhat.sim.ub.r1l0, labels = "-", cex = 1.5)
segments(x0 = 4, x1 = 4, y0 = yhat.sim.lb.r0l0, y1 = yhat.sim.ub.r0l0 )
text(x = 4 , y = yhat.sim.lb.r0l0, labels = "-", cex = 1.5)
text(x = 4 , y = yhat.sim.ub.r0l0, labels = "-", cex = 1.5)